Samenvatting

Om te laten zien dat GIS mogelijk is in R, gebruiken we wat open geodata van Den Haag, en maken een map.

Inleiding

Den Haag heeft een aantal kaarten openbaar gemaakt. Wij gebruiken de kaart van blokverwarming, energielabels (p5), en combineren deze met open satellietdata van NASA, en een open street map.

Data

Eerst downloaden en upzippen we de shapefiles van geoportaal-ddh.

# file download & unzip
download.file("http://geoportaal-ddh.opendata.arcgis.com/datasets/9cca53a6a4094a0b80e964b181a484ad_3.zip","blokverwarming.zip", mode="wb")
unzip("blokverwarming.zip")
download.file("http://geoportaal-ddh.opendata.arcgis.com/datasets/ed833f69fb82463ea63d37be00c356c3_0.zip", "energielabels.zip", mode="wb")
unzip("energielabels.zip")

# inputfileblokverwarming
f_blokverwarming <- "blokverwarming_DH.shp"
# inputfile enegielabels
f_energielabels <- "Energielabels_Postcode_5_niveau_Den_Haag_2016.shp"

Kaart Energielabels

We lezen de file in met readOGR van package rgdal. En daarna plotten we de polygons met wat kleurtjes, met behulp van base plot.

require(rgdal)
# de layers uit de shapefile
fl_energielabels <- ogrListLayers(f_energielabels) ## het is er slechts 1
# lees de layer in met readOGR
energielabels <- readOGR(f_energielabels,layer=fl_energielabels[1])
## OGR data source with driver: ESRI Shapefile 
## Source: "Energielabels_Postcode_5_niveau_Den_Haag_2016.shp", layer: "Energielabels_Postcode_5_niveau_Den_Haag_2016"
## with 844 features
## It has 10 fields
## Integer64 fields read as strings:  OBJECTID FREQUENCY labelCodeR MaxSimpTol
# plot de kaart
plot(energielabels,col=colorRampPalette(c("brown", "grey"))(5))

We willen natuurlijk graag laten zien wat de gemiddelde energielabels zijn per p5 gebied. Daarvoor gebruiken we de library tmap. Daarmee kun je vergelijkbaar met ggplot een kaart samenstellen.

require(tmap)
tm_shape(energielabels) +
  tm_polygons("MEAN_label", style="jenks", alpha=.5, palette=colorRampPalette(c("green", "red"))(5)) +
  tm_compass(type="arrow", position=c("right", "top"), fontsize = 2 ) + 
  tm_scale_bar()

## Kaart Blokverwarming Nu proberen we de file met blokverwarming in te lezen met readOGR.

# blokverwarming <- readOGR(f_blokverwarming)   # geeft een error

Die geeft een error, omdat blokverwarming een multiploint shapefile is. Dat formaat kan niet worden ingelezen met package gdal. We gebruiken st_read uit package sf. (NB readOGR heeft de voorkeur, mits die bruikbaar is natuurlijk.) Vervolgens plotten we de punten met behulp van base plot.

require(sf)
blokverwarming <- st_read(f_blokverwarming)
## Reading layer `blokverwarming_DH' from data source `C:\Users\Elisabeth\Desktop\Fourpoints\Coursera\GIS_in_R_voorbeeld_Den_Haag\blokverwarming_DH.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## geometry type:  MULTIPOINT
## dimension:      XY
## bbox:           xmin: 4.222274 ymin: 52.01797 xmax: 4.398967 ymax: 52.11275
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
plot(blokverwarming$geometry,pch='.')

## Kaart Energielabels met Blokverwarming Met behulp van base plot kunnen we de punten van blokverwarming over de p5 van energielabels heen.

# wat gevogel met x en y coordinaten uit het sf object halen
x_blok <- rep(1.1,length(blokverwarming$geometry[[1]])/2)
y_blok <- rep(1.1,length(blokverwarming$geometry[[1]])/2)
for (i in 1:length(blokverwarming$geometry[[1]])/2)
{
  x_blok[i] <-  blokverwarming$geometry[[1]][i,][1]
  y_blok[i] <-  blokverwarming$geometry[[1]][i,][2]
  i=i+1
}
blokverwarming_xy <- as.data.frame(cbind(x_blok,y_blok))
# plot
plot(energielabels,col=colorRampPalette(c("brown", "grey"))(5))
points(blokverwarming_xy,pch='.',col="blue")

We willen natuurlijk liever zo’n mooie tm_map plot maken.

blokverwarming_sp <- SpatialPoints(coords=blokverwarming_xy, proj4string = CRS(proj4string(energielabels) ) )
tm_shape(energielabels) +
  tm_polygons("MEAN_label", style="jenks", alpha=.5, palette=colorRampPalette(c("green", "red"))(5)) +
  tm_shape(blokverwarming_sp) +
  tm_dots(col="blue",legend.show = FALSE, legend.is.portrait = FALSE ) +
  tm_compass(type="arrow", position=c("right", "top"), fontsize = 2 ) + 
  tm_scale_bar()

## Open Street Map als referentie Alleen de vormen van p5 gebieden zijn eigenlijk niet genoeg voor een kaart. Om een basisreferentie van Den Haag/Nederland te hebben nemen we een OSM als achtergrond. We kunnen bijvoorbeeld ggmap gebruiken.

require(ggmap)
# Den Haag longitude and latitude
DH_loc = c(lon = 4.3, lat = 52.05 )
DH_basemap <- get_map(location = DH_loc, zoom = 12 )
ggmap(DH_basemap)

Maar met de instelling ‘tmap_mode(“view”)’ kun je een in-/uitzomobare kaart maken waar meteen al een basemap ondergeplot wordt. (Alleen kan dan het kompas niet laten zien worden.) Als je helemaal inzoomt kun je precies zien welk kadasterveld blokverwarming heeft.

tmap_mode("view")
tm_shape(energielabels) +
  tm_polygons("MEAN_label", style="jenks", alpha=.8, palette=colorRampPalette(c("green", "red"))(5)) +
  tm_shape(blokverwarming_sp) +
  tm_dots(col="blue",size=0.01) +
  tm_scale_bar()

Aantal panden met blokverwarming per p5

Een veel uitgevoerde bewerking bij het maken van kaarten is spatial join. Hierbij worden twee verschillende kaarten of layers samengevoegd. We voeren nu een spatial join uit van de p5 gebieden met energielabels, en de punten van de blokverwarming. En we maken een kaart met kleurindex op basis van aantal panden met blokverwarming.

a=rep(99999,length(energielabels@polygons))
for ( i in 1:length(energielabels@polygons))
{
  a[i]=sum(point.in.polygon(blokverwarming_xy$x_blok,
                   blokverwarming_xy$y,
                   energielabels@polygons[[i]]@Polygons[[1]]@coords[,1],
                   energielabels@polygons[[i]]@Polygons[[1]]@coords[,2]))
  i=i+1
}
energielabels$N_blokverwarming<-a
tm_shape(energielabels) +
  tm_polygons("N_blokverwarming", style="jenks", alpha=.8, palette=colorRampPalette(c("lightyellow", "brown"))(10)) +
  tm_shape(blokverwarming_sp) +
  tm_dots(col="blue",size=0.01) +
  tm_scale_bar()

Laten we ze naast elkaar plotten om ze te vergelijken.

tmap_mode("plot")
plot1 <- tm_shape(energielabels) +
  tm_polygons("MEAN_label", style="jenks", alpha=.8, palette=colorRampPalette(c("green", "red"))(5))+
  tm_compass(type="arrow", position=c("right", "top"), fontsize = 2 ) +
  tm_scale_bar()
plot2 <- tm_shape(energielabels) +
  tm_polygons("N_blokverwarming", style="jenks", alpha=.8, palette=colorRampPalette(c("lightyellow", "brown"))(10)) +
  tm_compass(type="arrow", position=c("right", "top"), fontsize = 2 ) +
  tm_scale_bar()
require(grid)
grid.newpage()
pushViewport(viewport(layout=grid.layout(1,2)))
print(plot1, vp=viewport(layout.pos.col = 1))
print(plot2, vp=viewport(layout.pos.col = 2))

Het lijkt niet gecorreleerd te zijn, maar we testen het lekker toch. Eerst kijken we of energielabel normaal verdeeld is. Helaas is dat niet zo, ook niet na een logtransformatie.

hist(energielabels$MEAN_label)

hist(log10(energielabels$MEAN_label))

## Raster Satelliet Data Via de NASA website heb ik rasterdata gedownload van de coordinaten van Den Haag. Deze data is niet direct te downloaden, maar de file staat in de github. Het is een .hdf file, hier is echt gdal voor nodig. Dit installeren op Windows had nog wat voeten in de aarde, maar het is uiteindelijk gelukt. Ik ben erg benieuwd hoe dat op andere besturingssystemen gaat.